ICM 2002 • Vol. Ill • 1-3 



Adaptive Finite Element Methods 
for Partial Differential Equations* 



R. Rannachert 

Abstract 

The numerical simulation of complex physical processes requires the use of 
economical discrete models. This lecture presents a general paradigm of deriv- 
ing a posteriori error estimates for the Galerkin finite element approximation 
of nonlinear problems. Employing duality techniques as used in optimal con- 
trol theory the error in the target quantities is estimated in terms of weighted 
'primal' and 'dual' residuals. On the basis of the resulting local error indica- 
tors economical meshes can be constructed which are tailored to the particu- 
lar goal of the computation. The performance of this Dual Weighted Residual 
Method is illustrated for a model situation in computational fluid mechanics: 
the computation of the drag of a body in a viscous flow, the drag minimization 
by boundary control and the investigation of the optimal solution's stability. 
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1. Introduction 

Suppose the goal of a simulation is the computation or optimization of a certain 
quantity J(u) from the solution u of a continuous model with accuracy TOL , by 
using the solution Uh of a discrete model of dimension N, 

A(u) = 0, A h (u h ) = 0. 

Then, the goal of adaptivity is the optimal use of computing resources, i.e., minimum 
work for prescribed accuracy, or maximum accuracy for prescribed work. In order 
to reach this goal, one uses a posteriori error estimates 

\J(u)-J(uh)\ « ri(uh) ■= ^ PK{u h )uj K , 

A'6T h 
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in terms of the local residuals pxiuh) of the computed solution and weights ujk 
obtained from the solution of a linearized dual problem. In the following, we will 
describe a general optimal control approach to such error estimates in Galerkin finite 
element methods. For earlier work on adaptivity, we refer to the survey articles jl()| . 
and 0. The contents of this paper is based on material from 0, [Hj and 0, 
where also references to other recent work can be found. 

2. Paradigm of a posteriori error analysis 

We develop a general approach to a posteriori error estimation for Galerkin 
approximations of variational problems. The setting uses as little assumptions as 
possible. Let X be some function space and L(-) a differentiable functional on X . 
We are looking for stationary points of L(-) determined by 

L\x){y) = VyeX, 

and their Galerkin approximation in finite dimensional subspaces Xh C X , 

L'{x h )(y h ) = \fy h € X h . 

For this situation, we have the following general result: 

Proposition 1 There holds the a posteriori error representation 

L(x) - L{x h ) = \L'{x h ){x-y h ) + R h , (2.1) 
for arbitrary yh € Xh- The remainder Rh is cubic in e :— x — Xh, 

Rh := \ \ L'" (x h + se)(e, e, e) s(s — 1) ds. 
Jo 

Proof We sketch the rather elementary proof. First, we note that 

L(x) - L(x h ) = / L'(xh+se)(e) ds 
Jo 

- \{L'{x h ){e) + L'(x)(e)} + \L'(x h ){e). 
Since Xh is a stationary point, 

L'(x h )(e) = L'(x h )(x-y h ) + L' '(x h )(y h -x h ) = L'(x h )(x-y h ), y h G X h . 
Finally, using the error representation of the trapezoidal rule, 

f(s) ds - i{/(0) + /(l)} = | f f'(s)s(s-l) ds, 

Jo 

completes the proof. Notice that the derivation of the error representation l|2.1|l does 
not assume the uniqueness of the stationary points. But the a priori assumption 
Xh — > x (h— »0) makes this result meaningful. 
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3. Variational equations 

We apply the result of Proposition ^ to the Galerkin approximation of varia- 
tional equations posed in some function space V, 

a{u){^) = Vip e V (3.1) 

Suppose that some functional output J(u) of the solution u is to be computed 
using a Galerkin approximation in finite dimensional subspaces Vh C V, 

a(u h ){ijj h ) = W/, € V h . (3.2) 

The goal is now to estimate the error J(u)—J(uh) ■ To this end, we employ a formal 
Euler-Lagrange approach to embed the present situation into the general framework 
laid out above. Introducing a 'dual' variable z ('Lagrangian multiplier'), we define 
the Lagrangian functional C(u,z) :— J(u) — a(u){z) . Then, stationary points 
{u, z} EVxV of £(■, •) are determined by the system 

rU \i i\ / J'{ u ){f) ~ d{v){<p,z) 1 n wr n 
£ {u,z)((p,ijj) = | _a(u){^j) J V{<P,ip}- 

The corresponding Galerkin approximation determines {uh, Zh} € VhxVh by 

Set x :— {u, z}, a;^ := {uh, Zh}, and L(a;) := C(u, z) . Then, 

J(u)-J(u h ) = L(x) + a(u)(z) - L(x h ) - a(u h )(z h )- 

Proposition 2 With the 'primal' and 'dual' residuals 

p(uh)(-) ■= -a(u h )(-), 
p*{z h ){-) := J'(u h )(-)-a'(u h )(;z h ), 

there holds the error identity 

J{u)-J{u h ) = \p{u h ){z-ijj h ) + \p*(z h )(u-ip h ) + Uh, (3.3) 

for arbitrary iph^h £ Vh ■ The remainder TZh is cubic in the primal and dual 
errors e u :— u — Uh and e z := z — Zh ■ 

The evaluation of the error identity 1|3.3[) requires guesses for primal and dual 
solutions u and z which are usually generated by post-processing from the ap- 
proximations Uh and Zh , respectively. The cubic remainder term TZh is neglected. 
We emphasize that the solution of the dual problem takes only a 'linear work unit' 
compared to the solution of the generally nonlinear primal problem. 
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4. Optimal control problems 

Next, we apply Proposition ^ to the approximation of optimal control prob- 
lems. Let V be the 'state space' and Q the 'control space' for the optimization 
problem 

J(u, q) -> min! a(u){ip) + b(q, ip) = Vip 6 V. (4.1) 

Its Galerkin approximation uses subspaces Vh X Q/i C V xQ as follows: 

J(u h , q h ) -> min! a(u h )(ij) h ) + b(q h , i>) = V?/^ £ V^. (4.2) 

For embedding this situation into our general framework, we again employ the Euler- 
Lagrange approach introducing the Lagrangian functional C(u, q, z) := J(u, q) — 
A(u)(z) — B(q, z) . Corresponding stationary points x := {it, q, z} e X := VxQxV 
are determined by the system ('first-order optimality condition') 

JL(u>q){<P) ~ a'(u)(<p,z) 1 

J' q (u,q)( X )-Hx,z) }=Q Vfe Xl # (4-3) 
-a(«)W-%,^) J 

The Galerkin approximation detrmines Xh '■= {uh, qh, z%\ 6 X% := Vh x x Vh in 
finite dimensional subspace V/, C F, Qf, C Q by 

J' u {uh-,qh){<Ph) - a'(uh)(<Ph,Zh) ] 

Jq(u h ,qh)(Xh) -b(xh,z h ) >=0 V{^ ft ,Xft, V^}- (4-4) 

-a(uh)(iph) -b(qh,iph) ) 

For estimating the accuracy in this discretization, we propose to use the natural 
'cost functional' of the optimization problem, i.e., to estimate the error in terms of 
the difference J(u, q)—J(uh, qh) ■ Then, from Proposition^ we immediately obtain 
the following result: 

Proposition 3 With the 'primal', 'dual' and 'control' residuals 

P*( z h)(-) ■= J' u ( u h,qh){-) - a'(u h )(-,z h ), 
P q ilh){-) := J' q (u h ,q h )(-)-b(; Zh ), 
P( u h)(-) ■= -a(u h )(-) -b(q h ,-), 

there holds the a posteriori error representation 

J(u,q)-J(u h ,q h ) = \p*{z h ){u-(ph) + \p q {qh){q-Xh) , . 

(4.5) 

+ 1 p(u h )(z-%l) h ) + TZh, 

for arbitrary tfh, iph € V% and \h € Qh- The remainder TZh is cubic in the errors 
e u := u-u h , e q := q-qu , e z := z-z h . 
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We note that error estimation in optimal control problems requires only the 
use of available information from the computed solution {uh, Ihi^h} , i- e -, n o extra 
dual problem has to be solve. This is typical for a situation where the discretiza- 
tion error is measured with respect to the 'generating' functional of the problem, 
i.e. the Lagrange functional in this case. In the practical solution process the 
mesh adaptation is nested with an outer Newton iteration leading to a successive 
'model enrichment'. The 'optimal' solution {u^ pt ,(?^ pt } obtained by the adapted 
discretization may satisfy the state equation only in a rather week sense. If more 
'admissibility' is required, we may solve just the state equation with an better dis- 
cretization (say on a finer mesh) using the computed optimal control g^ p * as data. 

5. Eigenvalue problems 

Finally, we apply Proposition ^ to the Galerkin approximation of eigenvalue 
problems. Consider in a (complex) function space V the generalized eigenvalue 
problem 

a(u,ip) = Xm(u,ip) \/tp GV, A g C, m(u,u) = 1, (5-1) 

where the form a(-, •) is linear but not necessarily symmetric, and the eigenvalue 
form 77i(-, •) is symmetric and positive scmi-dcfinit. The Galerkin approximation is 
defined in finite dimensional subspaces Vh C V , 

a(u hl ip h ) = X h m(u h , ip h ) Vtln £ 14, X h 6 C, m(u h , u h ) = 1. (5.2) 

We want to control the error in the eigenvalues X — Xh ■ To this end, we embed 
this situation into the general framework of variational equations by introducing 
the spaces V := VxC and Vh '•= V/,xC, consisting of elements U := {u, A} and 
Uh '■= {Uh, Xh} , and the semi-linear form 

A(U)(^) := Xm(u,ip)-a(u,ip) + ~p{m(u,u) - l}, * = {ip,n} £ V. 

Then, the eigenvalue problem (|5.1|l and its Galerkin approximation 1)5. 2J1 can be 
written in the compact form 

A(f7)(*) = VI eV, (5.3) 
A(U h )(*h) = V^eV,,. (5.4) 

The error in this approximation will be estimated with respect to the functional 

J($) := fim((p,ip), 

where J(U) — X since m(u, u) = 1 . The corresponding continuous and discrete 
dual solutions Z = {z,ir} € V and Zh = {zh,^h} £ V/i are determined by the 
problems 

A'(U)(9, Z) = J'(U)($) V$ e V, (5.5) 
A'(U h )($ h , Z h ) = J'(U h )($ h ) V$ h e V h . (5.6) 
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A straightforward calculation shows that these dual problems are equivalent to the 
adjoint eigenvalue problems associated to (|5.1|) and 15.21) . 

a(ip, z) — 7r m(ip, z) V<yj G V, m(u,z) = l, (5-7) 

a((p h , z h ) = ir h m((p h , z h ) V(p h e V h , m{u h , z h ) = 1. (5.8) 

Then, application of Proposition ^ yields the following result: 

Proposition 4 With the 'primal' and 'dual' residuals 

p(u h , X h )(-) ■= a(u h , •)-A /l 77i(u, l , •), 
p*{z hl ir h )(-) := a(-,Zh)-Trhm(-,z h ), 

there holds the a posteriori error representation 

X-X h = ^p(u h ,X h )(z-ip h ) + ±p*(z h ,ir h )(u-(p h ) -U h , (5.9) 

for arbitrary iph, fh S Vh, with the remainder term 

Tlh = \{X-\h)m{v-v h ,z-Zh)- 

We note that in Proposition 0] no assumption about the multiplicity of the 
approximated eigenvalue A has been made. In order to make the error represen- 
tation H5.9H meaningful, we have to use a priori information about the convergence 
{Aft, Vh} — ► {A, v} as h — > . The simultaneous solution of primal and dual eigen- 
value problems naturally occurs within an optimal multigrid solver of nonsymmetric 
eigenvalue problems. Further, error estimates with respect to functionals J(u) of 
eigenfunctions can be derived following the general paradigm. Finally, in solving 
stability eigenvalue problems A'(u)v = XA4v , we can include the perturbation of 
the operator A'(uh) ~ A'(u) in the a posteriori error estimate of the eigenvalues. 

6. Application in fluid flow simulation 

In order to illustrate the abstract theory developed so far, we present some 
results for the application of 'residual-driven' mesh adaptation for a model problem 
in computational fluid mechanics, namely 'channel flow around a cylinder' as shown 
in the figure below. The stationary Navier-Stokes system 

^„ ):= {-"A« + »-V« v +V P j =0 

determines the pair u := {v,p} of velocity vector v and scalar pressure p of a 
viscous incompressible fluid with viscosity v and normalized density p = 1 . The 
physical boundary conditions are v|rv iKid — , v\r ia — v m > and vd n v — np|r out = , 
i.e., the flow is driven by the prescribed parabolic inflow v ln . The Reynolds number 
is Re = ^— ^ = 20 , such that the flow is stationary. 
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Let the goal of the simulation be the accurate computation of the effective 
force in the main flow direction imposed on the cylinder, i.e. the so-called 'drag 
coefficient', 

J(u) := Cdrag = 1 in , 2n / n T {2vT-pI)eids 1 

where S is the surface of the cylinder, D its diameter, and r = ^{Vv+ Vv T ) the 
strain tensor. In practice, one uses a volume-oriented representation of Cdrag ■ 

Here, we cannot describe the standard variational formulation of the Navier- 
Stokes problem and its Galerkin finite element discretization in detail but rather 
refer to the literature; see 0,01 and the references therein. 

In the present situation the primal and dual residuals occuring in the a poste- 
riori error representation (|3.3() have the following explicit form: 

p{u h )(z-z h ) := ^2 [(R h ,z v -z v h ) K + {r h ,z v -z v h ) dK + (z p -zl,V-v h ) K + •■•}, 

K€T h 

p*(z h )(u-u h ) := ^2 {( R h v ~ v h)K + (ri,v-v h ) aK + (p-p h ,V-zl) K + . . .\, 
KeT h 

with the cell and edge residuals defined by 

/ + vAvh—Vh-Vvh- Vp, 

j + vAzl+v h -Vzl-VvZzl+V-v h zl-Vz p h , 

( \[vd n v h -np h ), if T(£dfl 1 
\ -vd n v h +nph, if T c r out , (= else) J ' 

f \[ud n zl+n-v h zl-z p h n], if T £ dfl \ 
\ -vd n z v h -n-v h zl + z p n, if T C T out , (= else) J ' 

where [. . . ] denots the jump across edges T , and '. . . ' stands for terms representing 
errors due to boundary and inflow approximation as well as stabilization. 

Practical mesh adaptation on the basis of the a posteriori error estimates pro- 
ceeds as follows: At first, the error functional may have to be regularized according 
to J{u) — J(u) + 0(TOL) . Then, after having computed the primal approximation 
Uh , the linear discrete dual problem is solved: 

(A'(uh)*z h ,iph) = J'{u h ){fh) V^eV^. (6.1) 

The error estimator is localized, 77^ = X)/xeT h ^ K > anc ^ approximation of the weights 
are computed by patch-wise higher-order interpolation: (z— Zh)\K ~ iJ-2h z h~ z h)\K ■ 



Rh\K '■— 

n h\K ■— . 

rh\r ■= 

* 

r h\r ■= 
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Finally, the current mesh is adapted by 'error balancing' rjx ~ f? u /#{-ftT 6 T^} . In 
the following, we show some results which have been obtained using mesh adaptation 
on the basis of the Dual Weighted Residual Method ('DWR method'). 

6.1. Drag computation (from |3j) 

The drag is computed on meshes generated by the DWR method and by an 'ad 
hoc' refinement criterion based on smoothness properties of the computed solution. 

Table 1: Results for drag computation on adapted meshes (1%-error in bold face). 



Computation of drag 


L 


N 


Qirag 


^drag 




4 


984 


5.66058 


l.le-1 


0.76 


5 


2244 


5.59431 


3.1e-2 


0.47 


6 


4368 


5.58980 


1.8e-2 


0.58 


6 


7680 


5.58507 


8.0e-3 


0.69 




oo 


5.57953 








Figure 1: Refined meshes by 'ad hoc' strategy (top) and DWR 
method (bottom) 

6.2. Drag minimization (from [4J) 

The drag coefficient is to be minimized by imposing a pressure drop at the 
two outlets r,; above and below the cylinder. In this case of 'boundary control' the 
control form is given by b(q 7 ip) —( ( lj r i-ip v )r 1 ur2 ■ 

Table 2: Uniform refinement versus adaptive refinement for Re = 40 . 



Uniform refinement 


Adaptive refinement 


N 


J drag 


N 


J drag 


10512 


3.31321 


1572 


3.28625 


41504 


3.21096 


4264 


3.16723 


164928 


3.11800 


11146 


3.11972 
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Figure 2: Velocity of the uncontrolled flow (top), controlled flow 
(middle), corresponding adapted mesh (bottom) 

6.3. Stability of optimized flows (from |8J) 

We want to investigate the stability of the optimized solution u opt = {v° pt , p opt } 
by linear stability theory. This is a crucial question since in the present case the 
optimal solution is obtained by a stationary Newton iteration which may converge 
to physically unstable solutions. In this context, we have to consider the non- 
symmetric eigenvalue problem for u := {v,p} £ V and A € C : 

ptx _/ -^Aw + w opt -Vw + wVw opt + Vp )_ J D 1 
*V )U -\ V-v J I J 

If the real parts of all eigenvalues are positive, Re A > , then the (stationary) 
base flow {v opt ,p° pt } is considered as stable (but with respect to possibly only 
very small perturbations) . We find that the optimal solution is at the edge of being 
unstable. 




Figure 3: Streamlines of real parts of the 'critical' eigenfunction 
shortly before the Hopf bifurcation and after, depending on the 
imposed pressure drop 
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Pressure drop Pressure drop 

Figure 4: Real and imaginary parts of the critical eigenvalue as 
function of the control variable 
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